TP_Estadistica_Espacial
Mora Assereto Farroni
2023-04-07
Librerías
library(geoR)## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.9-2 (built on 2022-08-09) is now loaded
## --------------------------------------------------------------library(geoR,lib.loc="PATH_TO_geoR")
library(sp)
library(sf)## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUElibrary(tmap)
#> Error : invalid version specification '1,5'
library(e1071)
library(spdep)## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(gstat)
library(caret)## Loading required package: ggplot2## Loading required package: latticelibrary(nlme)
library(ade4)##
## Attaching package: 'ade4'## The following object is masked from 'package:spdep':
##
## mstreeContexto
Este conjunto de datos contiene el contenido de magnesio en la capa de suelo de 20 a 40 cm en 178 ubicaciones dentro de una determinada área de estudio dividida en tres subáreas. La elevación en cada ubicación también se registró.La primera región suele estar inundada durante la temporada de lluvias y no se usa como área experimental. Los niveles de calcio representarían el contenido natural en la región. La segunda región ha recibido fertilizantes hace un tiempo y está típicamente ocupada por campos de arroz. La tercera región ha recibido fertilizantes recientemente y se usa con frecuencia como área experimental.
Formato
El set de datos original es un data frame con 178 observaciones y 10 variables, de las cuales para este trabajo solo son consideradas las siguientes:
east este-oeste (east-west) coordenadas en metros
north norte-sur (north-south) coordenadas en metros
mg2040 contenido de calcio en la capa de suelo de 20-40cm, medida en mmolc/dm3
Estos datos están disponibles en el paquete geoR (Paulo J. Ribeiro Jr and Peter J. Diggle (2016)).
Para realizar los análisis exploratorios a seguir, utilizaremos el paquete geoR e inicialmente necesitaremos crear un objeto del tipo geodata, como sigue.
Recoleción de datos
Importamos los datos para R
datos <- read.csv("~/Regresion Avanzada/Grupo10.csv", header=TRUE)
View(datos)# Mostrar de forma compacta la estructura de un objeto R arbitrario
str(datos) ## 'data.frame': 178 obs. of 4 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ east : int 5710 5727 5745 5764 5781 5799 5817 5837 5873 5907 ...
## $ north : int 4829 4875 4922 4969 5015 5062 5109 5161 5254 5347 ...
## $ mg2040: int 16 21 22 12 36 22 25 25 28 18 ...head(datos)## X east north mg2040
## 1 1 5710 4829 16
## 2 2 5727 4875 21
## 3 3 5745 4922 22
## 4 4 5764 4969 12
## 5 5 5781 5015 36
## 6 6 5799 5062 22class(datos)## [1] "data.frame"Transformamos el data set a un objeto de geodata
datos_geo <- as.geodata(datos, coords.col=2:3, data.col=4)class(datos_geo)## [1] "geodata"Para hacer el análisis de datos espaciales de manera más práctica, es común transformar el sistema de coordenadas geográficas en un sistema de coordenadas cartesianas, como el sistema de proyección UTM. Esto permite que las distancias entre los sitios se expresen en metros en lugar de grados, lo que facilita el análisis. Por lo tanto, el primer paso en el análisis de datos espaciales sería convertir las coordenadas geográficas en coordenadas cartesianas (UTM), en este caso no es necesario dado que las variables north y east que representan las coordenadas se encuentran expresadas en metros.
Análisis Exploratorio Espacial
En el análisis exploratorio de datos geoestadísticos es importante explorar la distribución de la variable. Esto se puede hacer mediante estadísticas descriptivas y gráficos de distribución de frecuencias. Si se supone que los datos tienen una distribución normal, estas medidas pueden ayudar a verificar los supuestos. Una distribución simétrica y similar a la de una variable normal se da cuando la media y la mediana son casi iguales y el coeficiente de asimetría es inferior a 1. La distribución de la variable también puede ayudar en la depuración de datos raros.
Breve resumen descriptivo de los datos
#Da un rápido resumen del objeto geodata cargado.
summary(datos_geo)## Number of data points: 178
##
## Coordinates summary
## east north
## min 4957 4829
## max 5961 5720
##
## Distance summary
## min max
## 43.01163 1138.11774
##
## Data summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.00000 21.00000 25.50000 25.58427 30.00000 43.00000Gráficamos los valores en las posiciones espaciales (distinguiendo según cuartiles), los datos frente a las coordenadas y un histograma de los datos. Los gráficos de dispersión de los datos frente a las coordenadas nos pueden ayudar a determinar si hay una tendencia
plot(datos_geo)
Generamos un gráfico con las posiciones de los datos (y por defecto con el tamaño de los puntos proporcional al valor):
points(datos_geo, col = "gray", pt.divide = "equal")
Transformación
datos_sf <- st_as_sf(datos, coords=c("east","north"),
crs=32721)head(datos_sf)## Simple feature collection with 6 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 5710 ymin: 4829 xmax: 5799 ymax: 5062
## Projected CRS: WGS 84 / UTM zone 21S
## X mg2040 geometry
## 1 1 16 POINT (5710 4829)
## 2 2 21 POINT (5727 4875)
## 3 3 22 POINT (5745 4922)
## 4 4 12 POINT (5764 4969)
## 5 5 36 POINT (5781 5015)
## 6 6 22 POINT (5799 5062)Eliminación de outliers e inliers
Las funciones hist() y boxplot() realizan
gráficos de histogramas y box-plots, respectivamente. Sus múltiples
argumentos permiten la edición de cada gráfico
par(mfrow = c(1, 2))
hist(
datos_sf$mg2040,
col = 'grey',
nclass = 20,
main = "Histograma",
ylab = 'Frecuencia Relativa',
xlab = 'Contenido Magnesio (mmolc/dm3)'
)
boxplot(
datos_sf$mg2040,
col = 'grey',
ylab = 'Contenido Magnesio (mmolc/dm3)',
main = "Box-Plot",
ylim = c(0, 50)
)
La función skewness() del paquete e1071
permite calcular el coeficiente de asimetría. Existen 3 fórmulas para su
cálculo (por defecto usa el tipo 3)
skewness(datos_sf$mg2040)## [1] 0.1099413En el histograma se observa que los datos tienen una distribución similar a una normal. No se ve una asimetría hacia algun lado en particular sino una disribución uniforme, también puede advertirse con los estadísticos con el coeficiente de asimetría el cual es de 0,109. En el gráfico box-plot no se observan valores extremos de la variable que se encuentren principalmente por encima de la media + 3SD
Las siguientes instrucciones calculan y crean objetos para la media, el DE y los límites superior (media + 3DE) e inferior (media - 3DE) con los que pueden detectarse los outliers y eliminar restos valores caso que sea necesario.
Media <- mean(datos_sf$mg2040)
DE <- sd(datos_sf$mg2040)
LI <- Media-3*DE
LS <- Media+3*DE# Rango de valores de outliers
(datos_sf$mg2040 > LS |
datos_sf$mg2040 < LI)## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [37] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [49] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [97] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [145] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [157] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [169] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSEComo no existe ningun valor dentro de este rango no es necesario la eliminación de outlieres
Detección de Inliers
La identificación y eliminación de inliers requiere de la
creación de una matriz de ponderación espacial. La función
dnearneigh() se utiliza para identificar el vecindario de
cada sitio. Para ello es necesario calcular la distancia espacial entre
los puntos para lo cual se usa la sintaxis $geom que
permite acceder a las coordenadas del objeto datos_sf. En
este ejemplo, se consideran datos vecinos a aquellos que se encuentran a
una distancia Euclídea de 0 a 100 m. La función
nb2listw() transforma el objeto vecindarios
que contiene las distancias a una matriz de pesos estandarizados por
filas (style = "W"). Este objeto es denominado
pesos_sp. Para generar la matriz de pesos espaciales es
necesario que todos los puntos tengan al menos un vecino, caso contario
la función nb2listw() genera un error advirtiendo este
hecho.
# 8345 es de hacer la distancia norte y este al cuadrado raiz seria el valor máximo
vecindarios <- dnearneigh(datos_sf$geom,
d1 = 0, d2 = 100)
# Muestra la cantidad de vecinos para la posicion
card(vecindarios)## [1] 3 5 7 6 6 7 6 6 5 6 7 5 5 6 8 8 9 9 9 7 7 7 8 8 9
## [26] 6 7 7 9 10 10 11 10 8 8 9 8 7 9 10 11 10 10 11 9 9 9 8 8 5
## [51] 6 7 10 11 9 10 9 8 8 8 8 5 5 7 10 12 11 10 10 10 9 7 7 6 7
## [76] 10 12 11 10 9 9 8 7 5 3 3 8 10 12 11 10 10 10 10 8 4 2 9 10 12
## [101] 11 11 10 10 9 9 7 5 6 8 10 10 10 9 9 6 7 5 8 11 11 11 10 9 10
## [126] 6 4 2 5 10 10 11 10 10 10 10 6 3 4 7 10 11 10 12 10 9 9 4 3 5
## [151] 5 6 10 11 11 10 9 9 8 6 1 5 2 7 10 10 10 8 7 5 1 3 3 6 8
## [176] 8 7 6# Hace un resumen de todos los valores
summary(vecindarios)## Neighbour list object:
## Number of regions: 178
## Number of nonzero links: 1416
## Percentage nonzero weights: 4.469133
## Average number of links: 7.955056
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9 10 11 12
## 2 3 7 4 15 17 21 22 25 41 16 5
## 2 least connected regions:
## 161 171 with 1 link
## 5 most connected regions:
## 66 77 89 100 144 with 12 linkspesos_sp <- nb2listw(vecindarios,
style = "W")La función localmoran() calcula el índice local de Moran
que permite identificar potenciales outliers espaciales. La
información referida al valor del índice local de Moran de cada punto se
encuentra en la columna Ii mientras que su significancia
estadística en la columna Pr(z < 0).
moranl <-
localmoran(datos_sf$mg2040,
pesos_sp,
alternative = "less")
head(moranl)## Ii E.Ii Var.Ii Z.Ii Pr(z < E(Ii))
## 1 0.86581506 -0.011967143 0.69358099 1.0539941 0.8540572
## 2 0.48460289 -0.002737869 0.09499218 1.5812073 0.9430847
## 3 0.39070018 -0.001673685 0.04103973 1.9368582 0.9736187
## 4 -1.06995511 -0.024040584 0.67628376 -1.2718378 0.1017154
## 5 -0.42053851 -0.014133570 0.40162630 -0.6412803 0.2606703
## 6 -0.08158934 -0.001673685 0.04103973 -0.3944842 0.3466118La identificación y eliminación de valores atípicos requiere la creación de una matriz de ponderación espacial. Para ello, se utiliza la función “dnearneigh()” para identificar el vecindario de cada sitio. Esto se logra calculando la distancia espacial entre los puntos, para lo cual se utiliza la sintaxis “$geom”, que permite acceder a las coordenadas del objeto “datos_sf”.
En este ejemplo, se consideran datos vecinos a aquellos que se encuentran a una distancia Euclidiana de entre 0 y 100 metros. La función “nb2listw()” transforma el objeto “vecindarios”, que contiene las distancias, en una matriz de pesos estandarizados por filas (style = “W”). Este objeto se denomina “pesos_sp”.
Es importante tener en cuenta que, para generar la matriz de pesos espaciales, es necesario que todos los puntos tengan al menos un vecino. En caso contrario, la función “nb2listw()” generará un error para advertir sobre este hecho.
moranp <-
moran.plot(
datos_sf$mg2040,
col = 3,
pesos_sp,
labels = F,
quiet = T,
xlab = "Contenido Magnesio",
ylab = "Contenido Magnesio Spatially Lagged"
)
Para visualizar en una tabla los puntos potencialmente influyentes y
sus estadísticos de diagnóstico se imprime el objeto
moranp. Datos con * en la columna inf se los considera como
influyente y por lo tanto posible outlier espacial.
summary(moranp)## x wx is_inf labels
## Min. :10.00 Min. :10.00 Mode :logical Length:178
## 1st Qu.:21.00 1st Qu.:23.38 FALSE:167 Class :character
## Median :25.50 Median :25.65 TRUE :11 Mode :character
## Mean :25.58 Mean :25.77
## 3rd Qu.:30.00 3rd Qu.:28.50
## Max. :43.00 Max. :32.30
## dfb.1_ dfb.x dffit cov.r
## Min. :-0.6807993 Min. :-0.3462534 Min. :-0.6992421 Min. :0.852
## 1st Qu.:-0.0229104 1st Qu.:-0.0242642 1st Qu.:-0.0474738 1st Qu.:1.006
## Median :-0.0022348 Median : 0.0011303 Median :-0.0012518 Median :1.016
## Mean :-0.0001416 Mean : 0.0001139 Mean : 0.0001092 Mean :1.011
## 3rd Qu.: 0.0247714 3rd Qu.: 0.0242995 3rd Qu.: 0.0733676 3rd Qu.:1.019
## Max. : 0.3771651 Max. : 0.6195409 Max. : 0.3847945 Max. :1.059
## cook.d hat
## Min. :5.000e-08 Min. :0.005640
## 1st Qu.:3.711e-04 1st Qu.:0.006374
## Median :2.449e-03 Median :0.008144
## Mean :6.454e-03 Mean :0.011236
## 3rd Qu.:6.341e-03 3rd Qu.:0.013070
## Max. :2.227e-01 Max. :0.044910Desde el objeto moranp puede extraerse una matriz de
valores lógicos (verdadero/falso) para los estadísticos diagnóstico que
identifican un punto como influyente o no.
influ <- moranp$is_inf
head(influ)## [1] FALSE FALSE FALSE TRUE FALSE FALSEEn la siguiente sentencia se adiciona al objeto datos_sf
los valores de los objetos moranl e influ, que
tienen información para detectar los outliers espaciales
detectados con el índice de Moran local y gráfico de Moran,
respectivamente.
datos_sf_1 <- cbind(datos_sf, moranl, influ)Posteriormente procedemos a eliminar los datos con Índice de Moran
Local negativo y estadísticamente significativos (p<0,05). La función
subset() selecciona datos que cumplen con cierta condición
lógica. El operador lógico or que indica que extraiga los datos
que cumplen con alguna de las dos condiciones. El nuevo objeto es
denominado datos_sf_2. Además, se crea un nuevo objeto que
tendrá los datos que han sido eliminados en este proceso
(inliers_ml).
datos_sf_2 <-
subset(datos_sf_1,
datos_sf_1[["Ii"]] >= 0 |
datos_sf_1[["Pr.z...E.Ii.."]] > 0.05)
inliers_ml <-
subset(datos_sf_1,
datos_sf_1[["Ii"]] < 0 &
datos_sf_1[["Pr.z...E.Ii.."]] < 0.05)summary(datos_sf_2$mg2040)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.00 21.00 25.50 25.56 30.00 43.00skewness(datos_sf_2$mg2040)## [1] 0.1121928par(mfrow = c(1, 2))
hist(
datos_sf_2$mg2040,
col = 'grey',
nclass = 20,
main = "Histograma",
ylab = 'Frecuencia Relativa',
xlab = 'Contenido Magnesio (mmolc/dm3)'
)
boxplot(
datos_sf_2$mg2040,
col = 'grey',
ylab = 'Contenido Magnesio (mmolc/dm3)',
main = "Box-Plot",
ylim = c(0, 50)
)
Luego de identificar y eliminar los inliers detectados con el índice de Moran y posteriormente con el gráfico de Moran, la nueva base de datos presenta 166 casos, es decir, se eliminaron 12 casos (6.7% de los datos) respecto a la base sin outliers.
Los estadísticos descriptivos de los datos depurados muestran una mejora en el coeficiente de asimetría (-0,19) lo cual se refleja en el histograma y box-plot.
Índice de Moran
Para la conformación de la matriz de ponderadores espaciales se
definieron los vecindarios de cada sitio mediante una red de conexión
construida en base a la distancia euclídea. Se consideraron sitios
vecinos a aquellos contiguos ubicados hasta 100 m de distancia. El
procedimiento es similar al empleado para el cálculo de índice de Moran
local. En este caso se agrega el argumento zero.policy=T
dentro de la función nb2listw() y moran.mc().
Esto permite que se genera la matriz de pesos espaciales sin la
restricción de que todos los puntos tengan al menos un dato vecino.
vecindarios <- dnearneigh(datos_sf_2$geom,
0, 100)
pesos_sp <- nb2listw(vecindarios,
style = "W",
zero.policy = TRUE)Para realizar el cálculo del Índice de Moran y determinar su
significancia estadística mediante simulación Monte Carlo, se utiliza
moran.mc(). Se debe especificar la variable en estudio, la
lista con los pesos de las ponderaciones espaciales y el número de
simulaciones.
i.moran <-
moran.mc(
datos_sf_2$mg2040,
listw = pesos_sp,
nsim = 999,
zero.policy = T
)
i.moran##
## Monte-Carlo simulation of Moran I
##
## data: datos_sf_2$mg2040
## weights: pesos_sp
## number of simulations + 1: 1000
##
## statistic = 0.37802, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greaterEstos resultados permiten concluir que existe autocorrelación espacial positiva (0,37802) y que esta es estadísticamente significativa (p=0,001).
Índice de Geary
geary.test(datos_sf_2$mg2040, listw = pesos_sp)##
## Geary C test under randomisation
##
## data: datos_sf_2$mg2040
## weights: pesos_sp
##
## Geary C statistic standard deviate = 8.6896, p-value < 2.2e-16
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.612803101 1.000000000 0.001985465El índice de Geary C es una medida de autocorrelación espacial que se utiliza comúnmente en análisis espaciales y de patrones espaciales. El resultado muestra que el valor de Geary C calculado para los datos es significativamente menor que el valor esperado (1.000000000), lo que indica una fuerte autocorrelación espacial en los datos. El valor de Geary C es 0.612803101, lo que sugiere que los valores similares están agrupados juntos en el espacio, mientras que los valores disímiles están más dispersos.
El estadístico desviación estándar de Geary C (8.6896) es alto y el valor p es extremadamente bajo (<2.2e-16), lo que indica que la probabilidad de obtener un valor tan bajo de Geary C por azar es extremadamente baja. Por lo tanto, se rechaza la hipótesis nula de que no hay autocorrelación espacial en los datos y se acepta la hipótesis alternativa de que existe una autocorrelación espacial significativa en los datos.
En resumen, el resultado del índice de Geary C sugiere que los valores similares están agrupados juntos en el espacio y que hay una fuerte autocorrelación espacial en los datos.
Análisis basado en semivariogramas opcion 1
Las semivariogramas empíricos pueden obtenerse usando la función
variogram()
Variograma nube
v <- variogram(mg2040~ 1,
datos_sf_2,
cloud=T,
cutoff=700)
plot(v)
Variograma
semiv_emp <- variogram( mg2040~ 1,
datos_sf_2,
cutoff = 600)
head(semiv_emp)## np dist gamma dir.hor dir.ver id
## 1 470 60.36988 26.33723 0 0 var1
## 2 656 108.13523 37.95655 0 0 var1
## 3 692 150.87729 43.01951 0 0 var1
## 4 575 181.80241 43.71304 0 0 var1
## 5 901 214.79382 44.33574 0 0 var1
## 6 1025 258.59180 46.49463 0 0 var1plot(semiv_emp,
main = "Contenido Magnesio (mmolc/dm3)",
xlab = "Distancia",
ylab = "Semivarianza")
En el caso anterior la fórmula utilizada (mg2040~1)
asume que el proceso es estacionario.
A continuación se ajusta un modelo de semivariograma teórico sobre el semivariograma empírico
Esférico
mod_esf <- fit.variogram(
semiv_emp,
vgm(0.6, "Sph", 200, 0.2))
mod_esf## model psill range
## 1 Nug 10.41296 0.0000
## 2 Sph 36.06021 196.9126plot(semiv_emp,
mod_esf,
main = "Contenido Magnesio (mmolc/dm3)",
xlab = "Distancia",
ylab = "Semivarianza")
El modelo que mejor ajusta será el de menor SCE. La función
attr() devuelve atributos de un objeto y puede usarse para
consultar la SCE del modelo ajustado.
attr(mod_esf, 'SSErr')## [1] 0.4988919Exponencial
mod_exp <- fit.variogram(
semiv_emp,
vgm(0.6, "Exp", 200, 0.2))
mod_exp## model psill range
## 1 Nug 0.00000 0.00000
## 2 Exp 48.63929 75.68158attr(mod_exp, 'SSErr')## [1] 0.2319651plot(semiv_emp,
mod_exp,
main = "Contenido Magnesio (mmolc/dm3)",
xlab = "Distancia",
ylab = "Semivarianza")
vgLine <- rbind(
cbind(
variogramLine(
mod_exp, maxdist = max(semiv_emp$dist)),
id = "Exponencial"),
cbind(
variogramLine(
mod_esf, maxdist = max(semiv_emp$dist)),
id = "Esférico")
)
ggplot(semiv_emp, aes(x = dist, y = gamma,
color = id)) +
geom_line(data = vgLine) +
geom_point() +
labs(
title = "Semivariograma experimental
y teórico ajustado") +
xlab("Distancia") +
ylab("Semivarianza") +
scale_color_discrete(
name = "Modelo",
labels = c("Esférico",
"Exponencial",
"Experimental"))
modelos <- fit.variogram(semiv_emp,
vgm(c("Exp", "Sph" ,"Gau" ,"Mat")))## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?modelos## model psill range kappa
## 1 Nug 0.00000 0.00000 0.0
## 2 Mat 48.63929 75.68156 0.5plot(semiv_emp,modelos)
attr(modelos, 'SSErr')## [1] 0.2319651Prediccion Espacial
La predicción espacial, es decir la predicción de valores de la variable en sitios del campo espacial donde no existen observaciones, usualmente se hace por el método kriging basándose en el semivariograma ajustado. Kriging proporciona el mejor estimador lineal insesgado del valor esperado para el sitio y un error de estimación conocido como varianza kriging. Esta varianza depende del modelo de semivariograma ajustado y de la ubicación en el espacio de los datos observados ya que son los datos observados en distintos sitios los que proveen información para aproximar el valor en el sitio sin dato.
limites <- read.csv("~/Regresion Avanzada/borders.csv", header=TRUE)
View(limites)
head(limites)## east north
## 1 5340 5800
## 2 5590 5690
## 3 5990 5690
## 4 5990 5100
## 5 5780 4800
## 6 5590 4800grid <- pred_grid(limites, by = 10)
grid <- polygrid(grid, bor = limites)
names(grid) <- c("x", "y")
grid <- st_as_sf(grid, coords = c("x", "y"),
crs = 32721)
plot(grid)
En resumen, estos pasos permiten generar una grilla regular de puntos dentro de un polígono determinado y transformarla en un objeto que puede ser utilizado en análisis de estadística espacial.
La función krige() del paquete gstat se utiliza para realizar interpolación kriging y simulaciones condicionales con diferentes métodos de predicción. En este caso, se utiliza la interpolación por kriging ordinario con el modelo de semivariograma matern estimado con geoestadística clásica. La función toma como argumentos la fórmula para especificar la estacionariedad del proceso, la base de datos, el objeto para la predicción, y la información del modelo de semivariograma teórico ajustado. Los argumentos opcionales nmin y nmax permiten realizar la interpolación en un contexto local.
# Kriging Ordinario
kriging_o <-
krige(
mg2040 ~ 1,
datos_sf_2,
st_as_sf(grid),
nmin = 7,
nmax = 25,
model = modelos
) ## [using ordinary kriging]#A continuación, se realiza la visualización de la predicción y su varianza.
predK_otm <- tm_shape(kriging_o) +
tm_dots("var1.pred", style = "cont",
title = "Predicción")
varK_otm <- tm_shape(kriging_o) +
tm_dots("var1.var", style = "cont",
title = "Varianza")
tmap_arrange(predK_otm, varK_otm)
Validación cruzada
La validación cruzada k-fold es un proceso que se utiliza para evaluar la precisión de un modelo de interpolación. En esta secciónn realizaremos los siguientes pasos:
Realizar la validación cruzada k-fold utilizando la función krige.cv() del paquete gstat.
Utilizar los argumentos fórmula (mg2040~1, lo cual especifica que es un proceso estacionario), la base de datos (datos_sf_2) y el modelo de semivariograma teórico ajustado (modelos) en la función krige.cv().
Definir el número de grupos (k) en los que se dividirá la base de datos utilizando el argumento nfold.
Fijar la semilla mediante la función set.seed() para obtener repetibilidad en los resultados.
Calcular estadísticos resumen como el error medio (ME), error cuadrático medio (MSE), media del cociente de la desviación cuadrática (mean squared deviation ratio, MSDR), raíz del error cuadrático medio (RMSE), la RMSE relativa a la media de los observados (RMSE_rel) y la correlación lineal entre los observados vs. predichos.
Utilizar estos estadísticos resumen para evaluar el desempeño del modelo de interpolación kriging.
Generar un gráfico de correlación entre los observados y predichos para visualizar la calidad del ajuste del modelo.
Realizar la validación cruzada, seleccion de k y escrbiri kriging
set.seed(17)
valcru <-
krige.cv(
mg2040 ~ 1,
datos_sf_2,
modelos,
nfold = 10,
nmin = 7,
nmax = 25
)Calcular los estadisticos de resumen y predicciones
ME <- mean(valcru$residual)
MSE <- mean(valcru$residual ^ 2)
MSDR <- mean(valcru$zscore ^ 2)
RMSE <- sqrt(mean(valcru$residual ^ 2))
RMSE_rel <-
sqrt(mean(valcru$residual ^ 2)) /
mean(valcru$observed) * 100
r <- cor(valcru$observed,
valcru$observed - valcru$residual)
tabla <- data.frame(ME, MSE, RMSE,
RMSE_rel, MSDR, r)
tabla## ME MSE RMSE RMSE_rel MSDR r
## 1 -0.2437268 27.86059 5.278313 20.65048 1.04758 0.6196965plot(
valcru$observed,
valcru$observed - valcru$residual,
xlab = "Observados",
ylab = "Predichos"
)